rm(list=ls())

setwd("~/Dropbox/COLOMBIA_PRIVATE/WRITING/Machine_Learning/aaaa-submitted/replication-archive/to-post/sims")

sim1.sc <- read.csv("sim1.csv")
sim5.sc <- read.csv("sim5.csv")
sim10.sc <- read.csv("sim10.csv")

# Results now are as ATT. Rescale to RIE based 
# on decomposition in paper
# Marginal probability of Z=1 is 0.31

sim1 <- sim1.sc*(-.31)
sim5 <- sim5.sc*(-.31)
sim10 <- sim10.sc*(-.31)

rmse1 <- sqrt(apply(apply(sim1, 
		1,
		function(x)(x[-1]-x[1])^2),
		1,
		mean, na.rm=T))
rmse5 <- sqrt(apply(apply(sim5, 
		1,
		function(x)(x[-1]-x[1])^2),
		1,
		mean, na.rm=T))
rmse10 <- sqrt(apply(apply(sim10, 
		1,
		function(x)(x[-1]-x[1])^2),
		1,
		mean, na.rm=T))

bias1 <- (apply(sim1, 2, mean, na.rm=T)-apply(sim1, 2, mean, na.rm=T)[1])[-1]
sd1 <- apply(sim1, 2, sd, na.rm=T)[-1]
bias5 <- (apply(sim5, 2, mean, na.rm=T)-apply(sim5, 2, mean, na.rm=T)[1])[-1]
sd5 <- apply(sim5, 2, sd, na.rm=T)[-1]
bias10 <- (apply(sim10, 2, mean, na.rm=T)-apply(sim10, 2, mean, na.rm=T)[1])[-1]
sd10 <- apply(sim10, 2, sd, na.rm=T)[-1]
sd.t <- sd(c(sim1[,1], sim5[,1], sim10[,1]), na.rm=T)

bOut <- rbind(bias1/sd.t, bias5/sd.t, bias10/sd.t)
sOut <- rbind(sd1/sd.t, sd5/sd.t, sd10/sd.t)
rOut <- rbind(rmse1/sd.t, rmse5/sd.t, rmse10/sd.t)

pdf(file="sim-out.pdf", width=6, height=9)
par(mfrow=c(3,1))
par(pty="m")
par(mar=c(5,5,2,1))
xMax <- 12
plot(	c(0,xMax),
		range(bOut),type="n",
		xlab="Number of noise covariates",
		ylab=bquote(Bias/sigma[psi]),
		axes=F)
axis(1, c(0,5,10))
axis(2, seq(-.2,.2, .05))
abline(h=0, lty="dashed")
points(c(0,5,10),
		bOut[,"OLS"], type="l")		
points(c(0,5,10),
		bOut[,"matching"], type="l")		
points(c(0,5,10),
		bOut[,"naiveIPW"], type="l")	
points(c(0,5,10),
		bOut[,"learnerIPW"], type="l")	
text(c(10,10,10,10),
	 bOut[3,],
	 c("OLS","Matching","Naive IPW","Ensemble IPW"),
	 pos=4)
box()

par(mar=c(5,5,2,1))
plot(	c(0,xMax),
		range(sOut),type="n",
		xlab="Number of noise covariates",
		ylab= bquote(S.E./sigma[psi]),
		axes=F,
		ylim=c(.85,1.1))
axis(1, c(0,5,10))
axis(2, c(.8, .9, 1, 1.1))
points(c(0,5,10),
		sOut[,"OLS"], type="l")		
points(c(0,5,10),
		sOut[,"matching"], type="l")		
points(c(0,5,10),
		sOut[,"naiveIPW"], type="l")	
points(c(0,5,10),
		sOut[,"learnerIPW"], type="l")	
text(c(10,10,10,10),
	 sOut[3,]+c(0,.005,-.005,0),
	 c("OLS","Matching","Naive IPW","Ensemble IPW"),
	 pos=4)
box()

par(mar=c(5,5,2,1))
plot(	c(0,xMax),
		range(rOut),type="n",
		xlab="Number of noise covariates",
		ylab= bquote(RMSE/sigma[psi]),
		axes=F,
		ylim=c(0,.25))
axis(1, c(0,5,10))
axis(2, seq(0,.25,.05))
abline(h=0, lty="dashed")
points(c(0,5,10),
		rOut[,"OLS"], type="l")		
points(c(0,5,10),
		rOut[,"matching"], type="l")		
points(c(0,5,10),
		rOut[,"naiveIPW"], type="l")	
points(c(0,5,10),
		rOut[,"learnerIPW"], type="l")	
text(c(10,10,10,10),
	 rOut[3,]+c(0,.005,-.005,0),
	 c("OLS","Matching","Naive IPW","Ensemble IPW"),
	 pos=4)
box()

dev.off()


